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Abstract 


The Standards of Fundamental Astronomy (SOFA) is a service provided by the International Astronomical Union 
that offers algorithms and software for astronomical calculations, which was released in two versions for FORTRAN 
77 and ANSI C, respectively. In this work, we implement the Python package PyMsOfa for SOFA service by three 
ways: (1) a Python wrapper package based on a foreign function library for Python (ctypes), (2) a Python wrapper 
package with the foreign function interface for Python calling C code (cffi) and (3) a Python package directly written 
in pure Python codes from SOFA subroutines. The package PyMsOfa has fully implemented 247 functions of the 
original SOFA routines released on 2023 October 11. In addition, PyMsOfa is also extensively examined, which is 
exactly consistent with those test examples given by the original SOFA. This Python package can be suitable to not 
only the astrometric detection of habitable planets from the Closeby Habitable Exoplanet Survey mission, but also 
for the frontier themes of black holes and dark matter related to astrometric calculations and other fields. The source 
codes are available via http://pypi.org/project/PyMsOfa/ and https: //github.com/CHES2023 /PyMsOfa. 
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1. Introduction 


The Standards of Fundamental Astronomy (SOFA)* service 
was established by the International Astronomical Union 
(IAU), which is an algorithm and program based on the 
fundamental astronomical models that aim to implement the 
TAU resolutions in an authoritative program. The SOFA service 
contains the implementation of the time system, the 
transformation of the coordinate system and the attitude of 
the Earth (e.g., the precession and nutation, closely related to 
International Earth Rotation Service), and astrometric para- 
meters. The routines were well written in strict compliance with 
IAU resolutions and have been updated as the resolutions were 
amended, with the latest version released on 2023 October 11. 

The SOFA service is officially released in two programming 
languages, i.e., FORTRAN 77 and ANSI C, whose advantages 
of high efficiency and fast running make them widely utilized in 
scientific computations. However, the Python language greatly 
benefits from its easy-to-learn, simple syntax, and tremendous 
built-in libraries in science and engineering, which motivates us 
to implement the full Python package PyMsOfa**® for SOFA, 
leading to direct, convenient and efficient routines being 
available. 
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There are several Python packages that contain SOFA 
features, such as Astropy,’ PyERFA,® and pysofa? (which 
contains only 186 subroutines from 247 SOFA service). 
However, in this work, the Python package is a complete 
implementation of all 247 functions of a recently released SOFA 
version in three ways via Python, i.e., wrapped with the ctypes 
library and the cffi interface, and written directly in pure Python 
codes based on SOFA subroutines. Thus the users of SOFA 
with FORTRAN 77 and ANSI C version can easily get started 
with this Python package, and achieve similar goals in the study. 

The paper is organized as follows: in Section 2 we briefly 
describe several major modules that SOFA service can achieve. 
Section 3 describes the features of the Python package and 
presents two examples of astrometric calculations. Finally, we 
make a concise summary in Section 4. 


2. Functions in the SOFA Service 


The routines of the SOFA service can be divided into four 
categories: basic calculation module, timescale and calendar 
module, coordinate system transformation module, and Earth 
attitude module. Each consists of dozens of routines. A brief 
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Table 1 
The Timescales Involved in SOFA Service 
Name Abb. Description 
International Atomic Time TAI Widely used time for the concrete realization of TT 
Coordinated Universal Time UTC The universal standard time 
Universal Time UTI Determined by the rotation of the Earth 
Terrestrial Time TT Coordinate time, can be used as the time argument of geocentric calendar 
former name: TDT 
Geocentric Coordinate Time TCG Coordinate time for the Geocentric Celestial Reference System (GCRS) 
Barycentric Coordinate Time TCB Coordinate time for the Barycentric Celestial Reference System (BCRS) 
Barycentric Dynamical Time TDB Time variable of the current solar system planetary calendar 


introduction of these routines and functions is provided below. 
More detailed descriptions can be found in the comments and 
SOFA documentation. 


2.1. Basic Calculation Module 


Such routines are very basic, mainly involving the 
processing of parameters, vectors and matrices. Its purpose is 
to conduct the calling between SOFA routines, so it does not 
fully include all vector and matrix operations. Similar 
functionality can be implemented in many other packages. 

These routines can realize the conversion of coordinates, 
such as the conversion between spherical coordinates to 
Cartesian ones (e.g., pymS2c and pymC2s). The routines 
about vectors perform the processing of velocity—position 
vectors that are required in other programs (e.g., pymTrxpv and 
pymPvu). For matrices, they contain functions such as rotation 
of a matrix (e.g., pymRx). In addition, the SOFA service also 
provides procedures for projective relations (e.g., pymTpxes), 
involving the conversion between spheres and planes. 


2.2. Timescale and Calendar Module 


The timescales involved in the SOFA service are shown in 
Table 1. Each timescale has a different use, and there is a slight 
time difference between them (Moyer 1981; Fairhead & 
Bretagnon 1990; Miiller & Jappel 2012). Some of the 
timescales have a linear relationship, while others have to be 
transformed according to the location (Seidelmann 1992; Soffel 
et al. 2003; McCarthy & Petit 2004). Figure 1 illustrates the 
functions involved in the conversion between timescales. In 
addition, the discontinuity of time caused by leap seconds also 
needs to be taken into account (e.g., pymDat) (Seidelmann 
1992). Therefore, when performing astronomical scientific 
calculations, it is necessary to choose an appropriate timescale. 

The Julian Date is an astronomical method that measures 
time over a long time span, which is widely adopted in the 
definition of a calendar. In the SOFA service, the Julian Date is 
described by two parameters, for a specific calendar of Julian 
Day JD = 2450123.7, the date can be denoted as 
JD1 = 2400000.5 and JD2 = 50123.2, respectively. 


2.3. Coordinate System Transformation Module 


The SOFA service include the conversion of different 
coordinate systems, such as geodetic coordinate, geocentric 
coordinate, horizon coordinate, hour angle system of 
coordinate, ICRS, ecliptic coordinate and galactic coordinate. 
The package provides a complete system for the astrometric 
parameters, ASTROM, which is composed of 30 components. 
These components include the time, solar system barycenter 
(SSB) to the observer, direction and distance relative to the 
Sun, barycentric observer velocity, reciprocal of Lorenz factor, 
bias-precession-nutation matrix, Terrestrial Intermediate Origin 
(TIO) locator (McCarthy & Petit 2004), polar motion, geodetic 
latitude, diurnal aberration (Klioner 2003), Earth rotation angle 
and refraction constants A and B (Crane 1976; Riieger 2002). 
Detailed calculation procedures are also provided for each 
component. The final ASTROM parameters are the astrometric 
parameters of different stations, independent of stars. Table 2 
lists several major basic coordinate systems involved in SOFA 
service, while Table 3 displays the functions that produce the 
components of ASTROM and the functions used to transform 
the stellar parameters. 

SOFA service also deals with the transitions between FK4, 
FK5 and the Hipparcos catalog (e.g., pymFk425, pymFk524 
and pymFkShip) (Aoki et al. 1983; Yallop et al. 1989; Mignard 
& Froeschlé 2000). These coordinate systems cover the 
fundamental coordinate systems in classical astrometry. In 
astronomical calculation, the change of coordinate system is 
inseparable from the change of position of the observation 
target to various observation stations (Green 1962). The SOFA 
service here gives a set of authoritative algorithms. 

In addition to this, the SOFA service still provides the 
transformation of star positions between FK4, FK5 and the 
Hipparcos catalog. 


2.4. Earth Attitude Module 


The SOFA service supplies an authoritative definition of the 
parameters involved in the Earth attitude and provides detailed 
algorithms. These parameters include the parameters of the 
objects in the solar system (e.g., the mean longitude of the 
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Figure 1. The functions involved in the conversion between the timescales in 
Table 1. 


Table 2 
The Fundamental Coordinate Systems Involved in SOFA Service 


Name Abb. Description 


International Celestial ICRS An ideal coordinate system that does not 


move with 


Reference System respect to distant extragalactic objects. 


Celestial Intermediate CIRS Geocentric reference system modified by 


precession 


Reference System and nutation based on GCRS. 


Geocentric Celestial GCRS Geocentric spacetime coordinate system 


within the 


Reference System framework of general relativity. 


Barycentric Celestial BCRS A set of spacetime coordinate systems 


of the 


Reference System solar system centroid. 


planets and their approximate heliocentric positions and 
velocities, e.g., pymPlan94) (Simon et al. 1994; Klioner 
2003), bias-precession-nutation matrix (Lieske 1979; Mathews 
et al. 2002; Wallace & Capitaine 2006), Earth rotation angle 
(Capitaine et al. 2000), celestial intermediate pole (CIP), 
celestial intermediate origin (CIO) (Capitaine et al. 2003a, 
2003b; Capitaine & Wallace 2006), etc. (e.g., pymPn00 and 
pymNum06a). There are some differences in these parameters 
under different models, so SOFA also gives the parameters 
under the IAU 1976 model (Lieske et al. 1977), IAU 1980 
model (Seidelmann 1982), IAU 2000 A&B model (Soffel et al. 
2003) and IAU 2006 model (Seidelmann et al. 2007). The 
calculations of the routines are also consistent with those in the 
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IAU resolution (Wallace & Capitaine 2006). Thus one can 
choose the required models and parameters with respect to 
different astronomical calculations. 


3. Python Feature of the Package 


For simplicity, here this Python package of the SOFA 
service is described based on the ctypes library, which is a 
foreign function library for Python that provides C compatible 
data types, and allows calling functions in the shared 
libraries.'° We integrated all 247 functions in a C file from 
the SOFA package, which is named sofa_a.c. A dynamic 
shared library file (libsofa_c.so) was compiled from three 
source files (sofa.h, sofam.h and sofa_a.c). By calling the 
shared library file and the ctypes library, we establish a Python 
interface for each SOFA function under ANSI C, thus SOFA 
service can be employed directly in Python programs. The 
codes are accessible at https://github.com/CHES2023/ 
PyMsOfa and http: //pypi.org/project /PyMsOfa/. 

This Python package differs from the ANSI C version in that 
the parameters involved in its functions vary in addition to the 
language used. In the original ANSI C version by SOFA, the 
input and output parameters of the function are included in the 
parameters of the calling function. In the Python version, the 
parameters of the function refer only to those input parameters. 
In addition, the parameters involved in each function, along with 
their types, units, etc., are briefly given in their comments. It can 
be easily queried when using each function by the package in 
the Python environment. In the following, we will show two 
applications for the usage of the Python package for SOFA. 

For the Closeby Habitable Exoplanet Survey (CHES) (Ji 
et al. 2022), a number of functions in the SOFA service are 
employed to simulate the observed images. Here we present 
two examples (see Figures 2 and 3). The first example 
illustrates the time expression in SOFA service as mentioned 
previously in that the time is usually expressed by two 
parameters. The given two examples can accurately calculate 
the coordinates of the target star and reference stars based on 
the Gaia Data Release 3 (DR3), respectively, at the time of 
observation (on the plane). From the projection theorem, the 
coordinates in the focal plane are given to calculate the major 
observed quantity in the CHES mission, i.e., the angular 
distance between them in the focal plane. Such functions can 
further provide references to relevant parameters in the CHES 
mission (Ji et al. 2022; Tan et al. 2022). Figure 4 shows a 
simulation image with the above-mentioned two functions 
(pymPmsafe, pymTpxes), which is in good agreement with the 
produced image from ESASky. 

In addition to the detection of exoplanets, this package is 
also available for various astrometric observations, such as the 
detection of black holes, dark matter, etc. Figure 5 reproduces 


10 The cffi interface version and pure Python codes are also available via the 
GitHub website. 
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import PyMsOfa 
ral, dec1=0. 08788796644747092, -1. 132188403551008 
pmrl, pmd1=8. 27454225597756e-06, 5. 647882757242327e-06 
px1=116. 18257854098105e-3 
rv1=9. 283571 
epla, ep1b=2400000. 5, 57388. 450625000056 
ep2a, ep2b=2400000. 5, 60991. 538933193 
ra2, dec2, pmr2, pmd2, px2, rv2=PyMsOfa. pymPmsafe (ral, decl, pmrl, pmdl, pxl, 
rvl, epla, eplb, ep2a, ep2b) 
print (ra2, dec2) 
0. 08796958189264041 -1. 1321326881050646 
Figure 2. pymPmsafe: Example for the calculation of the R.A. and decl. at the given epoch. 
import PyMsOfa 
a, b=0. 0902535601897608, —1. 1308652476489707 
a0, b0=0. 08796958189264041, —1. 1321326881050646 
xi, eta=PyMsOfa. pymTpxes(a, b, a0, b0) 
print (xi, eta) 
0. 0009726944659440729 0. 001266436096867058 
Figure 3. pymTpxes: Example for the calculation of rectangular coordinates on the focal plane. 
Table 3 
The Transformation between Basic Coordinate Systems 
Observer Transformation Parameters Produce* Application” 
space ICRS > GCRS pymApcs pymApcs13 pymAtciq* pymAticq* 
geocentric ICRS «+ GCRS pymApcg pymApcg13 
terrestrial ICRS + CIRS pymApci pymApcil3 
ICRS +> Observed pymApco pymApcol13 pymAtciq* pymAticq* 
pymAtiog pymAtoiq 
CIRS « Observed pymApio pymApio13 pymAtioq pymAtoiq 
Notes. 


è The number 13 at the end of the functions indicates the use of internal functions to compute various ephemerides. 
P The symbol (*) indicates that the function has multiple variations, such as pymAtciqn and pymAtciqz. 


the prediction of the motion of the star orbiting a black hole carefully processed to achieve an accurate orbital motion under 
using PyMsOfa. The orbital parameters in the simulations are the gravitation from a black hole (Figure 5). This will 
adopted from El-Badry et al. (2023). Using this Python contribute to a more detailed investigation and analysis of the 
package, the proper motion and parallax of the star can be physical properties of black holes. 
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Figure 4. Left panel: The simulation image of the target star (* zet Tuc) in the CHES mission. The positions of the stars are based on the data from Gaia DR3 and 
converted to the predetermined observation time (JD2460992.04) and observation position (L2) by the relevant function in PyMsOfa. Right panel: The image of the 


same field of view obtained from ESASky. 
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Figure 5. The predicted motion of Gaia BH1’s photocenter on the sky during 5 yr. Left panel: The total motion due to proper motion, parallax and orbital motion. The 
dash-dotted line isolates the contribution due to parallax and proper motion alone. Right panel: The predicted orbital motion, with parallax and proper motion removed. 


4. Summary 

In this work we describe the functionality of the Python 
package PyMsOfa for SOFA service in three ways: (1) a Python 
wrapper package based on ctypes library, (2) a Python wrapper 
package with the cffi interface and (3) a Python package directly 
written in pure codes from SOFA algorithms. They implement 
SOFA service in a Python interface in various ways but each 
package can have the same usage. The purpose of this work is to 
enable a wide variety of usages of the authoritative algorithms in 
the SOFA service with the Python interface. The Python 
package has fulfilled the necessary requirements all 247 
functions of SOFA service written in ANSI C, containing basic 
calculation, timescale and calendar, coordinate system 


transformation, Earth attitude and astrometric parameters. This 
Python package is thoroughly and correctly tested with those 
original examples given by SOFA, which runs stably on Linux, 
macOS and Windows operating systems. 

Astrometry is an ancient branch of astronomy, which 
primarily aims at the investigation of the positions of celestial 
bodies. With quick development, astrometry can supply a clear 
understanding from the solar system to the entire universe on 
the basis of ground-based and space-based observations. 
Nowadays, with improvement in the accuracy of astrometry, 
besides direct measurements of stellar positions and velocities, 
indirect measurements can be extensively conducted to provide 
vital implications for the frontier themes of habitable planets, 
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black holes, dark matter, etc. The SOFA service integrates the 
basic theory of astrometry and is a convenient tool to conduct 
astrometric related calculations and can be applied to diverse 
research fields. This Python package was originally intended to 
facilitate the use of SOFA service for various astronomical 
calculations. Thanks to the excellent, comprehensive features 
of SOFA service, this package can also be further used for 
related astronomical calculations in Python. 
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